home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Cream of the Crop 20
/
Cream of the Crop 20 (Terry Blount) (1996).iso
/
math
/
alged34.zip
/
ALGEDSRC.ZIP
/
ALGEXPR.C
< prev
next >
Wrap
C/C++ Source or Header
|
1996-06-06
|
20KB
|
787 lines
/*--------------------------------------------------------------------
Alged: Algebra Editor henckel@vnet.ibm.com
Copyright (c) 1994 John Henckel
Permission to use, copy, modify, distribute and sell this software
and its documentation for any purpose is hereby granted without fee,
provided that the above copyright notice appear in all copies.
*/
#include "alged.h"
/*--------------------------------------------------------------------
associate
rotate the association on add,sub,mul,div
a+b-c => a-c+b => b+a-c => b-c+a
note, we make use of the fact that ADD+1 = SUB and ADD is even.
*/
void associate(node *p) {
node *a,*b;
int opr;
opr = p->kind;
b = p->rt;
if (opr==DIV && b->kind==MUL) { /* special handle for bisected */
p->rt = b->rt;
b->rt = b->lf;
b->lf = p->lf;
p->lf = b;
b->kind = DIV;
b = p->rt;
}
if (opr==ADD || opr==MUL || opr==SUB || opr==DIV) {
a = p;
while ((a->lf->kind|1) == (opr|1)) {
a->kind = a->lf->kind;
a->rt = a->lf->rt;
a = a->lf;
}
a->kind = opr;
if (opr&1) { /* subtr or divide */
a->rt = b;
}
else {
a->rt = a->lf;
a->lf = b;
}
}
else if (opr==EQU) { /* commute equality */
p->rt = p->lf;
p->lf = b;
}
}
/*--------------------------------------------------------------------
commute
commute on add,sub,mul,div
e.g. a/b ==> b^(-1)/a^(-1)
*/
void commuteNOTUSED(node *p) {
node *a,*b;
a = p->lf;
b = p->rt;
if (p->kind==MUL || p->kind==ADD || p->kind==EQU) {
p->lf=b; p->rt=a;
}
else if (p->kind==DIV || p->kind==SUB) { /* non-commutative */
if (a->kind==p->kind+1 &&
b->kind==p->kind+1 &&
a->rt->kind==NUM &&
b->rt->kind==NUM) { /* has coefficients */
p->lf=b; p->rt=a;
a->rt->value = -a->rt->value;
b->rt->value = -b->rt->value;
}
else { /* need to make coefficients */
p->lf = newoper(p->kind+1);
p->rt = newoper(p->kind+1);
p->lf->lf = b;
p->lf->rt = newnum(-1);
p->rt->lf = a;
p->rt->rt = newnum(-1);
}
}
}
/*--------------------------------------------------------------------
a^x * a^y ==> a^(x+y)
x^a * y^a ==> (xy)^a
*/
int expjoin(node *p) {
int i,r=0;
node *a,*b;
for (i=0; i<p->nump; ++i)
r+=expjoin(p->parm[i]);
a = p->lf;
b = p->rt;
if (p->kind==MUL || p->kind==DIV) {
if (a->kind==EXP && b->kind==EXP &&
equal(a->lf,b->lf)) { /* a^x * a^y = a^(x+y) */
freetree(a->lf);
a = a->rt;
freenode(p->lf);
b->kind = p->kind-2;
p->kind = EXP;
p->lf = b->lf;
b->lf = a;
}
else if (b->kind==EXP &&
equal(a,b->lf)) { /* a * a^y = a^(1+y) */
freetree(a);
a = newnum(1);
b->kind = p->kind-2;
p->kind = EXP;
p->lf = b->lf;
b->lf = a;
}
else if (a->kind==EXP &&
equal(a->lf,b)) { /* a^x * a = a^(x+1) */
freetree(b);
b = newnum(1);
a->kind = p->kind-2;
p->kind = EXP;
p->lf = a->lf;
a->lf = a->rt;
a->rt = b;
p->rt = a;
}
else if (a->kind==EXP && b->kind==EXP &&
equal(a->rt,b->rt)) { /* x^a * y^a = (xy)^a */
freetree(b->rt);
p->rt = a->rt;
a->rt = b->lf;
freenode(b);
a->kind = p->kind;
p->kind = EXP;
}
else return r;
++r;
}
return r;
}
/*--------------------------------------------------------------------
remove sub and div
x/y = x*y^-1
x-y = x+y*-1
*/
int nosubdiv(node *p) {
int i,r=0;
node *a,*b;
for (i=0; i<p->nump; ++i)
r+=nosubdiv(p->parm[i]);
a = p->lf;
b = p->rt;
if (p->kind==SUB || p->kind==DIV) {
++r;
--p->kind;
if (b->kind==NUM && p->kind==ADD) /* a-2 = a+(-2) */
b->value = -b->value;
else if (b->kind==NUM && !whole(b->value))
b->value = 1.0/b->value;
else {
a = newoper(p->kind+2); /* MUL or EXP */
a->lf = b;
a->rt = newnum(-1);
p->rt = a;
}
}
return r;
}
/*--------------------------------------------------------------------
bisect node
This converts expressions to canonical form in which
1. + * are left assoc, ^ is right assoc.
2. scope of / is increased, e.g. (a/b)*c ==> ac/b,
3. x+y*-2 is changed to x-y*2
3' x+(-2) is changed to x-2
4. x*y^-2 is changed to x/y^2
4' x*(0.5) is changed to x/2
5. scope of ^ is reduced
The name "bisect" means that each MUL clag is broken into two pieces the
numerator elements and the denominator ones.
*/
int bisect(node *p) {
int i,r=0;
node *a,*b;
for (i=0; i<p->nump; ++i)
r+=bisect(p->parm[i]);
a = p->lf;
b = p->rt;
switch (p->kind) {
case ADD:
case MUL:
/*--------------------------------------------------------------------
Add or Multiply
*/
if (b->kind==p->kind || /* a+(b+c) = a+b+c */
b->kind==p->kind+1) { /* a+(b-c) = a+b-c */
swingb;
i = b->kind;
b->kind = p->kind;
p->kind = i;
}
else if (p->kind==ADD && /* a+(-2) = a-2 */
b->kind==NUM &&
b->value < 0) {
p->kind = SUB;
b->value = -b->value;
}
else if (b->kind==p->kind+2 && /* a+b*-2 = a-b*2 */
b->rt->kind==NUM &&
b->rt->value < 0) {
++p->kind;
b->rt->value = -b->rt->value;
}
/*--------------------------------------------------------------------
Add only
*/
else if (p->kind==ADD &&
aop(b->kind)) { /* a+(b-c) = (a+b)-c */
swingb;
p->kind = b->kind;
b->kind = ADD;
}
/*--------------------------------------------------------------------
Multiply only
*/
else if (p->kind==MUL &&
a->kind==p->kind+1) { /* a-b+c = a+c-b */
p->rt = a->rt;
a->rt = b;
++p->kind;
--a->kind;
}
else if (p->kind==MUL &&
a->kind==p->kind+2 && /* b*-2+a = a-b*2 */
a->rt->kind==NUM &&
a->rt->value < 0) {
++p->kind;
p->lf = b;
p->rt = a;
a->rt->value = -a->rt->value;
}
else break;
++r; break;
case SUB:
case DIV:
/*--------------------------------------------------------------------
Subtract or divide
*/
if (b->kind==p->kind+1 && /* a-b*-2 = a+b*2 NOT NECES. */
b->rt->kind==NUM &&
b->rt->value < 0) {
--p->kind;
b->rt->value = -b->rt->value;
}
/*--------------------------------------------------------------------
Subtract only
*/
else if (p->kind==SUB && /* a-(-2) = a+2 */
b->kind==NUM &&
b->value < 0) {
p->kind = ADD;
b->value = -b->value;
}
else if (p->kind==SUB &&
aop(b->kind)) { /* a-(b+c) = (a-b)-c */
swingb;
p->kind = (ADD+SUB) - b->kind;
b->kind = SUB;
}
/*--------------------------------------------------------------------
Divide only
*/
else if (p->kind==DIV &&
b->kind==p->kind) { /* a-(b-c) = a+c-b */
p->lf = b;
p->rt = b->lf;
b->lf = a;
--b->kind;
}
else if (p->kind==DIV &&
a->kind==p->kind) { /* a-b-c = a-(b+c) */
swinga;
--a->kind;
}
else break;
++r; break;
case EXP:
/*--------------------------------------------------------------------
exponent
*/
if (a->kind==EXP) { /* (x^y)^z = x^(y*z) */
swinga;
a->kind = MUL;
}
else if (b->kind==NUM && /* a^(-2) = 1/a^2 */
b->value < 0) {
p->kind = DIV;
p->lf = newnum(1);
p->r